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ABSTRACT 

Nonlinear evolution is sometimes modeled by assuming there is a deterministic map- 
ping from initial to final values of the locally smoothed overdensity. However, if an 
underdense region is embedded in a denser one, then it is possible that its evolution is 
determined by its surroundings, so the mapping between initial and final overdensities 
is not as 'local' as one might have assumed. If this source of nonlocality is not accounted 
for, then it appears as stochasticity in the mapping between initial and final densi- 
ties. Perturbation theory methods ignore this 'cloud-in-cloud' effect, whereas methods 
based on the excursion set approach do account for it; as a result, one may expect 
the two approaches to provide different estimates of the shape of the nonlinear counts 
in cells distribution. We show that, on scales where the rms fluctuation is small, this 
source of nonlocality has only a small effect, so the predictions of the two approaches 
differ only on the small scales on which perturbation theory is no longer expected 
to be valid anyway. We illustrate our results by comparing the predictions of these 
approaches when the initial-final mapping is given by the spherical collapse model. 
Both are in reasonably good agreement with measurements in numerical simulations 
on scales where the rms fluctuation is of order unity or smaller. 

If the deterministic mapping from initial conditions to final density depends on 
quantities other than the initial density, then this will also manifest as stochasticity in 
the mapping from initial density to final. For example, the Zeldovich approximation 
and the Ellipsoidal Collapse model both assume that the initial shear field plays an 
important role in determining the evolution. We compare the predictions of these 
approximations with simulations, both before and after accounting for the 'cloud-in- 
cloud' effect. Our analysis accounts approximately for the fact that the shape of a cell 
at the present time is different from its initial shape; ignoring this makes a noticable 
difference on scales where the rms fiuctuation in a cell is of order unity or larger. 

On scales where the rms fiuctuation is 2 or less, methods based on the spherical 
model are sufficiently accurate to permit a rather accurate reconstruction of the shape 
of the initial distribution from the nonlinear one. This can be used as the basis for a 
method for constraining the statistical properties of the initial fiuctuation field from 
the present day field, under the hypothesis that the evolution was purely gravitational. 
We illustrate by showing how the highly non-Gaussian nonlinear density field in a 
numerical simulation can be transformed to provide an accurate estimate of the initial 
Gaussian distribution from which it evolved. 

Key words: methods: analytical - dark matter - large scale structure of the universe 



The present work is primarily concerned with the prob- 
ability distribution function (hereafter pdf; some authors 



1 INTRODUCTION 



prefer to call this the probability density function), which 
describes the probability that a randomly placed cell of 
specified shape and volume contains a certain amount of 
mass. In the best studied cosmological models, the pdf 
has a Gaussian form initially, but becomes increasingly 
positively skewed at later times. There are two methods 
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for estimating the evolution of the dark matter pdf. 
One is based on perturbation theory (jBernardeau 19941 



Protogeros fc Scherrer 1997 
Scherrer & Gaztanaga 2001 



'Gaztanaga & Cr oft 1999[ 
iBernardeau et al. 20d2| and 
the other is based on excursion set methods l|Sheth 1998p . 
The perturbation theory based calculation is not expected 
to be accurate on scales where the variance is of order 
unity or larger; in hierarchical models, this means that 
perturbation theory is not expected to be accurate on 
small scales. Hyperextended Perturbation Theory (HEPT) 
l|Scoccimarro fc Frieman 1999|) is expected to be valid in 
the larger variance (small cell) regime where standard 
perturbation theory breaks down; it provides explicit 
expressions for the moments of the pdf rather than a closed 
form expression. For the special case of clustering from 
white-noise Gaussian initial conditions, the excursion set 
method predicts exactly the same pdf as does HEPT, 
despite the fact that the methods used by the two ap- 
proaches are very different. Motivated by this coincidence, 
the purpose of the present note is twofold: first, to show 
that, for more general initial conditions, the excursion set 
approach actually makes rather similar predictions to those 
of perturbation theory on scales where the variance is small; 
the second is to understand why. 

In the analysis which follows, we distinguish between 
two steps in the calculation of the pdf. The first is the ap- 
proximation for nonlinear evolution which we will call the 
dynamics. The second is how this approximation is used to 
translate from the initial pdf to an evolved one, which we 
call the statistics. In the first half of this paper, we study 
approximations for the dynamics which are based on the as- 
sumption of spherical symmetry. In this case, the perturba- 
tion theory method provides what is, in effect, a monotonic, 
deterministic mapping between the initial and final overden- 
sities. Because the final overdensity at a specified position 
in space is determined solely by the initial value at that po- 
sition, this is sometimes also called a 'local' mapping, since 
values of the initial fluctuation field at other positions are 
assumed to not affect the mapping. 

The excursion set approach accounts for the fact that 
the evolution of a given region may actually be determined 
by less local surroundings. For example, consider the evolu- 
tion of an underdense region which is surrounded by a dense 
shell. If the shell is sufficiently dense, then it will eventu- 
ally collapse, crushing the smaller region within it. The lo- 
cal approximation would have predicted expansion rather 
than collapse for the smaller underdense region. Sheth & 
van de Weygaert (2004) call this the void-in-cloud problem, 
although it is clear that this is an extreme example of a more 
general 'cloud-in-cloud' problem. Clearly, in such cases, the 
mapping between initial and final overdensities is not as 'lo- 
cal' as perturbation theory assumes, and accounting for this 
'cloud-in-cloud' problem is likely to be more important for 
small 'clouds'. If not accounted for, this efi'ect will mani- 
fest both as stochasticity (since the same initial overdensity 
may map to many different final densities depending on the 
surroundings) and, perhaps, as a bias. The excursion set ap- 
proach provides an algorithm which accounts for this source 
of non-locality; it assumes that, once the correct large scale 
has been chosen, the mapping is deterministic. 

However, there is another source of non-locality which 
spherical evolution models ignore. This source plays a cru- 



cial role in models which account for the infiuence of the 
external shear field on the evolution. The simplest of these 
more sophisticated models for the dynamics is the Zeldovich 
approximation (Zeldovich 1970|) . Here, the nonlinear density 
is determined not just by the initial overdensity, but by two 
other quantities which are related to the surrounding shear 
field. These quantities also play an important role in exten- 
sions of the Zeldovich approximation (Maklcr ct al. 200lJ 
as well as in ellipsoidal collapse models (White fc S ilk 1979] 
[Bond fc Myers 1996[ ). In all these models, the nonlinear den- 
sity is a deterministic function of three quantities associated 
with the initial fluctuation fleld. In the context of perturba- 
tion theory models for the pdf, the mapping from initial 
density to final density will appear to be stochastic if the 
influence of the two other variables is not accounted for. In 
the excursion set approach, this stochasticity is in addition 
to that which derives from the cloud-in-cloud problem which 
is now associated with all three variables. 

To explore the consequences of this additional source 
of stochasticity we show the result of inserting the Ellip- 
soidal collapse model into the perturbation theory and ex- 
cursion set calculations. We do this in two steps: by showing 
the predictions when expanding the dynamics to first, and 
then to second order. To first order, the Ellipsoidal collapse 
model reduces to the Zeldovich approximation; hence, our 
perturbation theory discussion of the first order approxi- 
mation extends previous work on the Zeldovich approxima- 
tion HKofman et al. 19941 [Hui, Kofman fc Shandarin 2000| ). 
However, our analysis also accounts for another subtlety as- 
sociated with non-spherical collapse models — that the fi- 
nal shape of a patch is different from its initial shape 
(|Betancort-Rijo 1991 IPadmanabhan fc Subramanian 19931 



|Betancort-Rijo fc Lopez-Corredoira 2002p . We then present 
the results of the simplest excursion set treatment of this 
problem which accounts for the evolution in the volumes 
but not the shapes of cells. 

The perturbation theory and excursion set predictions 
associated with spherical dynamics are presented in Sec- 
tion (2] This section also contains a discussion of the Log- 
normal model. It includes a comparison of these predictions 
with measurements in simulations. The perturbation and ex- 
cursion set treatments of the Zeldovich and ellipsoidal col- 
lapse models are in Section [S] A final section summarizes 
our findings, and includes a discussion of the fact that local 
deterministic mapping models of nonlinear evolution moti- 
vate a simple method for reconstructing the shape of the 
initial pdf from that at late times. We illustrate the method 
using the spherical evolution model; we also show that, for 
the present purposes, the spherical model is a good enough 
approximation to the second order ellipsoidal collapse dy- 



2 DETERMINISTIC MAPPINGS FROM 
INITIAL TO FINAL DENSITY 

This section describes the perturbation theory and excur- 
sion set models of the pdf when the mapping from linear to 
nonlinear density is deterministic. In both, the variance of 
the initial density fiuctuation field when smoothed on scale 
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Rm plays an important role. It is denoted by 
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where Pi^ik) is the power spectrum of the initial field, W{x) 
is the Fourier transform of the smoothing window, and 
Rm = (3M/47rp)^/^ We wiU also use the quantity 



dlngj 
' dlnM' 



(2) 



For P(fc) oc fc", 7 = (n + 3). 

2.1 Perturbation theory-based methods 

These methods make three important assumptions. First, an 
initial overdensity can be mapped to an evolved density 
p. Because p at any given position in space is determined 
from 5l at the same position (in Lagrangian coordinate), 
the mapping is said to be 'local'. Because p depends on (5l 
and nothing else, the mapping is said to be 'deterministic'. 
Notice that these 'local' 'deterministic' assumptions make no 
mention of the scale on which they apply. In what follows, 
we will write the evolved density in a cell of volume V , in 
units of the background density, as 

M 

= i + (3) 

we hope that this slight abuse of notation will not cause 
confusion. Since the evolved density is clearly smoothed on 
scale V , the question is; on what smoothing scale should the 
initial overdensity 5l associated with p be defined? This is 
where the third key assumption is made: the appropriate 
scale is that which initially contains mass M. As a result of 
this assumption, the cumulative distributions of the evolved 
(Eulerian) pdf at fixed V and the initial (Lagrangian) pdf at 
fixed mass scale M are related as follows (see Section 5.4.3 
in Bernardeau et al. 2002): 



r°° M 



dx 



4l(M|V)/<tl(M) 



exp(-a;V2) 



(4) 

where M = dM p{M\V) M, ctl is given by equation JlJ, 
and the right hand side assumes that the initial distribution 
was Gaussian. Differentiating with respect to M yields 



exp[-(^L/crL)V2] d((?L/crL) 

dM 



M 



Since M/M = p, the expression above implies that 



p p(pjV) = exp 



Si dln(5L/<7L) 



2nal 



d In p 



(5) 



(6) 



The requirement that the left hand side of equation (j4| de- 
creases monotonically with increasing M (since we want 
p{M\V) > for all M) means that (5l/(Tl must increase 
monotonically with increasing Ad. As a result, the range 
of allowed ShiMlV) relations is constrained by the relation 
(jl(M), i.e., by the initial power spectrum. 

2.2 Excursion set method 



The excursion set model (|Sheth 199811 exploits the fact that 
a'l{M) is a monotonic function of A4, so the local collapse 



model Si^{M\V) defines a curve in the space of Si, versus 
al. In what follows, we wiU set S{M) = crl{M). Then, the 
excursion set model for the distribution of M in cells of size 
V is 

d\nS 



pMp\v) = SfiS\V) 



dlnp 



(7) 



where /(5'|V')dS' denotes the probability that a random 
walk with uncorrelated steps first crosses a barrier of height 
B{S\V) on scale S (where S denotes the variance in walk 
heights). The collapse dynamics is included in this solution 
of the statistical problem by setting B{S\V) = (5l(o"l)- It is 
by relating the counts in cells distribution to the first cross- 
ing distribution that the excursion set model accounts for 
the 'cloud-in-cloud' problem discussed in the Introduction. 
Note that this method allows a larger range of S-l{M\V) 
relations than does perturbation theory. 

We have used two methods for computing f{S\V)dS: 
one is a Monte-Carlo approach, where we simulate a large 
ensemble of random walks and count the distribution of first 
crossings, and the other is the analytical approximation of 
Sheth & Tormen (2002). This approximation sets 

B{S\Vf' 



.f(S\V)dS = \T{S\V)\ exp 



2S 



dS/S 



(8) 



where T{S\V) denotes the first few terms of the Taylor ex- 
pansion of B about S. For the barrier shapes of interest in 
this paper, including only the first two terms is sufficient, so 

dlnB{S\Vy 



T{S\V) ^ B{S\V) 



1 



dlnS 



(9) 



where the derivative is evaluated at S. Thus, 



p p{p\V) = exp 



exp 



B(S'|V)= 



25 
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1 - 
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(10) 



where the final expression uses the fact that B^/S = 
((Jl/ctl) . Comparison with perturbation theory (equa- 
tion |6]) shows that the only difference is in the Jacobian 
like term, which differs by a factor of ctl. 

2.3 Normalization 

Direct implementation of the two methods described above 
produces distributions which are guaranteed to have the 
correct mean value, J dpp{p\V) p = 1, but, in general, 
/ dpp{p\V) 7^ 1. (Strictly speaking, the analytic approxi- 
mation to the random walk does not integrate to unity; the 
Monte-Carlo solution, of course, does; nevertheless the ana- 
lytical formula gives a very good approximation.) To ensure 
that this integral also equals unity, one must define 

/oo 
dpp{p), p' = Np, and p'^p(p) = p'^p{p). (11) 



The quantities p' and p'{p') now satisfy both normalisation 
requirements: / dp'p'{p') = 1 and J dp'p'{p')p' — 1. 
Note that although this procedure is standard 
()Kofman et al. 19941 |Proto geros & Scherrer 19971 
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IQhta et al. 2004"|> . this procedure has no physical mo- 
tivation, other than to insure correct normahzation. 

The last of the expressions above shows that plots of 
p{p) versus log(p) can be turned into properly normalized 
plots simply by shifting along the log(p) direction. This is 
useful because, as equations ([6]) and ((TJ show, it is p{p) 
which is the fundamental prediction of both models for the 
statistics. We will use this in what follows. 



2.4 Example: The Lognormal distribution 

The standard form of the Lognormal distribution is 



2 exp[-(lnp + ^)V2cr2l P 

P PKP) = 7= 

V27r o- 



(12) 



where p — (t^/2 and = In(p^). This can be scaled to a 
Gaussian distribution in the variable (5l with mean zero and 
variance ctl by setting p= 1 + 5 oc exp(5L) so ln(l + 5) oc 5l 
where the constant of proportionality is set by requiring that 

(P) = 1- 

However, if one views the transformation p oc exp(5L) 
as a local deterministic mapping between the initial over- 
density 5l and evolved density p (see Coles & Jones 1991 
and Smith, Scoccimarro & Sheth 2007 for motivation), then 
equation (|12[) is not the form which one obtains from the 
perturbation theory-like argument we have just described. 
In this context, equation (|12|) corresponds to the Lagrangian 
rather than Eulerian space pdf. This curious fact does not 
appear to have received attention before. The distinction is 
important, because, for most power spectra of current inter- 
est, 5l(M|1/)/(jl(M) oc \n{M/V)/ai,{M) is not monotonic 
in M, so the method of Section [2.11 breaks. 



2.5 Example: Spherical collapse dynamics 

The spherical collapse model relates (5l and p; it is well ap- 
proximated by 



1 + 5^ 



5l 



(13) 



(Bernardeau 1994). The exact value of 5c depends on the 
background cosmology. It is 5c ~ 1.686 in an Einstein de- 
Sitter model, and 5c ~ 1.66 for the simulation results we 
present later. As a result, setting 5c = 5/3 is an excellent 
approximation which we will use to provide simple analytic 
examples of our results. Note that the model studied by 
Sheth (1998) was, in effect, 5c — 1; for white-noise initial 
conditions (initial spectral slope n = 0) , the resulting excur- 
sion set prediction is in exact agreement with HEPT. As we 
discuss in more detail later, this expression with Ssc = 3 is 
what the Zeldovich approximation yields for the collapse of 
a spherical perturbation. 

The perturbation theory prediction for the pdf associ- 
ated with the spherical collapse mapping (equation 1 13|) is 
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Figure 1. Comparison of equation IIA5I I (solid, short-dashed, dot- 
short dashed) and equation IIASI l (dotted, long-dashed, dot-long 
dashed) for av = 0.5, 1.0, 2.0 respectively. 



The corresponding excursion set expression is 

2 . 1 



pp{p\v) 



exp 



5l{p) 



Slip) 



2^lip) 



+ lSLip) 



(15) 



Note that this differs from the perturbation theory expres- 
sion only because it has a factor of 7/3 rather than 7/6. 
For |(5l(p)| ^ 1, or for 7 ^ 6/5bc ~ 18/5, this difference is 
negligible, and the two models will give very similar results. 
Demonstrating this agreement between the two methods is 
one of the central results of this paper. 

To illustrate the similarity, and to get a feel for the 
magnitude of the differences. Figure [1] shows the nonlinear 
pdf predicted by these models when 5c — 5/3 and the initial 
power spectrum was Pik) oc k~^^^. The three sets of curves 
are for av = 1/2, l,and 2 (narrowest to broadest distribu- 
tions). 



2.6 Comparison with simulations 

We compare our predictions with measurements made in 
the Very Large Simulation (VLS) (JV;gshida ot al. 2001). 
The simulation box represents a cube 479/i~^ Mpc on 
a side in a cosmological model where (flm, ^a, cs) = 
(0.3,0.7,0.7,0.9). It contains 512^ particles, so the associ- 
ated particle mass is 6.86 x lO^^/i"^ M©. 

We initially counted particles in cubic cells of side 
(479/37) /i~^Mpc; the volume of each cell is equivalent to 
that of a sphere of radius 8.03 ft~'^Mpc. In the results which 
follow, counts in larger cells were obtained by summing up 
counts in neighbouring cells. For cells near the boundary of 
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Figure 2. Comparison of tiie measured p{p) (filled squares) witli various splierical collapse based predictions for this quantity. Top: 
dot-dashed (green) curves show the standard Lognormal (equation [12}. Short-dashed (orange) and dotted (red) curves show perturbation 
theory before (equation [T4]l and after the normalization discussed in Section l2.3l (the dotted curve is simply shifted slightly to the right); 
solid (cyan) curves show the exact Monte-Carlo solution of the excursion set model, and long-dashed (blue) curves show the associated 
analytic approximation (equation I15I I. The nonlinear rms fluctuations of different models are shown in the following order: the VLS 
simulation (0.968 and 0.513), the Lognormal (1.069 and 0.560), the perturbation theory before normalization (0.921 and 0.523), the 
perturbation theory after normalization (1.004 and 0.539), the exact Monte-Carlo soultion of the excursion set model (1.134 and 0.559), 
and the analytic approximation (1.094 and 0.549); Bottom: the normalized distributions and the Lognormal divided by the measurements. 



the box, we used the fact that the simulation was run using 
periodic boundary conditions. 

Because one does not simulate the nonlinear evolution 
of a continuous density field, rather, one simulates the mo- 
tion of particles, the relation between the discrete particle 
pdf one measures in simulations is non-trivial. Typically, one 
assumes that the distribution of particle counts in cells is 

p{N\V)= j dMp{M\V)p{N\M), (16) 

where p{N\M) denotes the probability that a mass M 



is represented by A*' particles. If m denotes the particle 
mass, then the Poisson model assumes that p{N\M) = 
(M/m)^ exp{—A4/m)/N\. Our comparisons with the theo- 
retical models of the previous section assume that discrete- 
ness effects are negligible, so N/N — M/M. 

The top panel of Figure [2] shows the pdf of the dark mat- 
ter in cells of two different volumes (equivalent to spheres of 
radii 8 and 16 /i~^Mpc), chosen to have rms fluctuation val- 
ues of unity or less. Black squares show the measurements in 
the VLS, where p = 1 + Snl = N/N . Curves show the pre- 
dictions of the different models: the Lognormal is dot-dashed 
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Figure 3. Same as previous figure, but for smaller cells, in which the rms fluctuation is significantly larger than unity. The results from 
the exact Monte-Carlo solution of the excusrion set model are not included as the associated analytical approximation gives very similar 
results. 



(this model has one free parameter, the nonlinear variance, 
which we compute from the fitting formula in Smith et al. 
(2003), rather than using the value measured in the simu- 
lations); raw and normalized perturbation theory model are 
short-dashed and dotted (note that the normalised curve is 
simply shifted to the right of the unnormalised curves); and 
the excursion set model is solid (long-dashed curve shows 
the analytic approximation). The nonlinear variances a as- 
sociated with these various models are also shown (from left 
to right, top to bottom: VLS simulation, lognormal, unnor- 
malised perturbation theory, normalised perturbation the- 
ory, excursion set from Monte-Carlo, and excursion set from 
the approximation formula). 

To reduce the dynamic range, the bottom panel shows 
the predictions and Lognormal all ratioed to the measure- 
ments. This shows clearly that both the normalised pertur- 



bation theory and excursion set models are in better agree- 
ment with the simulations than is the Lognormal. We note in 
passing that Betancort-Rijo & Lopez- Corredoira (2002) pro- 
vide a relatively simple analytic expression for the nonlinear 
pdf. Their expression is, essentially, yet another prediction 
for the spherical evolution based pdf. It fares substantially 
worse than either of our models - it overpredicts the high 
density tail - so we do not show it here. 

Figure |3] shows that these spherical collapse based pre- 
dictions remain accurate even on scales where the rms fluc- 
tuation is significantly larger than unity. This is well beyond 
the regime where standard perturbation theory is expected 
to be valid, and indeed, the perturbation and excursion set 
predictions for the high density tails differ significantly. Per- 
turbation theory provides a substantially better description 
of the simulations when the rms fluctuation is of order 2 (left 
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panel) but is worse when the rms is larger (right panel). On 
the other hand, the average density within the virial radius 
of a dark matter halo is about 200 times the background 
density; thus, in the regime where In p > 5, the mapping of 
equation (|13[l is suspect. 

Figures [2] and |3] show that the spherical model produces 
rather accurate predictions for the shape of the nonlinear 
pdf, at least on scales where the nonlinear rms fluctuation is 
smaller than about 2. We show later that this can be used 
to motivate an algorithm which uses the nonlinear density 
field to provide an estimate of the initial pdf. 



3 STOCHASTIC MAPPINGS FROM INITIAL 
TO FINAL DENSITY 

In this section we replace the assumption of a spherical col- 
lapse in favor of the ellipsoidal collapse model. This evolu- 
tion is considerably more complex, and so we only present 
results to first and second order in the dynamics. We use the 
formulation of this model in which it reduces, to lowest or- 
der, to the Zeldovich Approximation ( |Bond fc Myers 1996[ ). 



3.1 The Zeldovich Approximation 

In the Zeldovich Approximation, the nonlinear density is a 
deterministic function of three numbers in the initial dis- 
tribution. These are the eigenvalues Xi of the deformation 
tensor, a 3 x 3 symmetric matrix whose elements are the 
second derivatives of the initial potential field. 

The initial density (5l is the sum of the three eigenvalues, 
whereas the nonlinear density is 



P = 



(17) 



(note that p ^ 1 + Ai when A <C 1). All three eigenval- 
ues are the same for a sphere, in which case A = (5l/3. In 
this case, the relation above reduces to that of the spheri- 
cal model from the previous section with 5c — 3. Similarly, 
5c = 2 is associated with regions where the smallest eigen- 
value is zero, and the other two are each equal to (5l/2. If 
one calls such an object a filament, then a sheet has two 
eigenvalues equal to zero and the third equal to 5l, so the 
evolution is given by equation (I13|) with 5c = 1. 
The evolved pdf of this model has 



pp{p\V)dp = 



(18) 
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Hui, Kofman & Shandarin 2000; 



Betancort-Rijo fc Lopez-Corredoira 2002[ ). 

For Gaussian initial conditions, the joint distribution 
of p(Ai, A2, Aslcr) is known (|Doroshkevich 1970|l . and it is 
straightforward to evaluate the integral above by Monte 
Carlo methods. In practice, this distribution is a function 
of Xi/a, so it is useful to think of the expression above as 



ppip\V) dp = / dAp(A|l) & P = - 



(19) 



Notice that if cr = (7l(pV) = av, then 



3 

p{p\V)dp^ J d\p{\\l) 11(1 



crvXi) — 1. 



(20) 



This choice is made by Hui et al. (2000); in this case, 
the normalization problems of Section 12.31 do not arise. 
However, the analysis of the previous sections suggests 
that this ignores the effects of smoothing, and that set- 
ting a = ai,{M) is almost certainly a better choice 



(|Betancort-Rijo 1991 IPadmanabhan fc Subramanian 19931 
|Betancort-Rijo fc Lopez-Corredoira 2002p . 

Notice that if any one of the eigenvalues equals unity 
then the density diverges, and that the density becomes neg- 
ative if one and only one or all three eigenvalues exceed 
unity. This happens often when cr > 1 (with probability 11% 
when (T = 1), thus restricting the use of this approximation 
to large scales where a is small. Accounting for the effects 
of smoothing will mitigate this somewhat, since ctm ^ cv 
when p ^ 1, where we have defined am = o-l{M) and 
(Jv = crL(pV). For simplicity, on the few occasions when 
the density does go negative, we take the absolute value of 
the rhs of equation (|17|) : there is no physical motivation for 
this choice, but, in practice, this happens sufficiently rarely 
that it does not affect our results. 

The expression above is still not entirely self-consistent. 
Although it tries to account for the fact that the final 
size of a region is different from its initial size, it does 
not account for the fact that the shape has also changed 
(|Betancort-Rijo 199 1 ; Padmanabhan fc Subramanian 19931 
|Betancort-Rijo fc Lopez-Corredoira 2002[ ). If the final vol- 
ume V is spherical, the initial volume was not; the analy- 
sis above does not account for this. However, this can be 
done in a relatively straightforward way by simply setting 
(jl(A), where, for Eulerian spheres V of radius -Re we require 
that _Re = Ri{i — Xi). This means that the delta function 
picks out ellipsoids in the initial field which contained mass 
M oc ni=i ^ I ^-nd which are now spheres of volume V which 
contain the same mass M. The function (tl(A) is the rms 
fiuctuation in the initial field when smoothed with a triaxial 
filter of shape given by {Ri, R2, R3); we have checked that 
it is well approximated by 



C^L (A) 



Sph 



(M) cxp 



B 



E 



In 



1- A, 
1 - A, 



(21) 



where B = 0.0486, and the subscripts ell and 
sph stand for ellipsoidal and spherical respectively 
( |Betancort-Rijo fc Lopez-Corredoira 2002| ). This, then, is 
our perturbation theory based estimate associated with the 
Zeldovich Approximation. 

Figure |3] shows how the shape depedence of the rms 
fiuctuation changes the predicted dark matter pdf. In the 
top panel the dotted curves are the predictions of pertur- 
bation theory using equation (|19|l with ai^{M) whereas the 
dashed curves use o"l(A). The difference between these two 
is significant, meaning the shape dependence is important, 
in the high density tail (ln(p) > 1) of the 8/i~^Mpc cells. 
The bottom panel compares the ratios of different models 
to the measurements. In neither cases does the Zeldovich 
Approximation agree well with the simulations: the differ- 
ence is more than a factor of two in low density regions 
(ln(p) < -2 in 8/i"^Mpc and ln(p) < -1.5 in 16/i"^ Mpc). 
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-2 2 -2 2 

ln(p) ln(p) 

Figure 4. Comparison of the measured p(p) (filled squares) and various Zeldovich approximation-based predictions and second 
order ellipsoidal model using perturbation-theory for this quantity. In the top panel red (dotted) and cyan (solid) curves show the 
normalized perturbation-theory and excursion-set predictions using the spherical linear variance approximation. Magenta (dash) shows 
the perturbation-theory prediction using the ellipsoidal variance approximation (equation [21} . Blue (dotted-dash) shows the normalized 
perturbation-theory predictions using second order ellipsoidal collapse model (equation [SOj . The numbers shown are the rms values 
for (from top to bottom) VLS simulation, excursion-set (spherical variance), perturbation-theory (spherical variance), perturbation- 
theory (Zeldovich approximation using ellipsoidal variance), and perturbation-theory (second ellipsoid collapse using ellipsoidal variance). 
Bottom panel shows the corresponding distributions divided by the measurements. 



It is rather more complicated to account approximately 
for the eflect of the change in shape in the excursion set 
approach. Our Monte Carlo algorithm works as follows. We 
generate (Ai,A2,A3) in each step of the walk following the 
procedure described by Sheth & Tormen (2002). The vari- 
ance associated with step n is S^"\ This has an associ- 
ated scale We are interested in patches which today 
are spheres of radius Re- (The extension to ellipsoids at 
the present time is straightforward). These are those initial 
patches which satisfy R^"\l — A^"') = Re, where A^"' is 
the value of Xi after n steps. We then want the largest mass 



associated with the various Rl" . Since mass is proportional 
to Yli^ij want the largest value of -R^"' for each i; in 
effect, we want the first crossing values Ui for each of the 
three barriers R^/^\l — A^"') = Re- Let /(ni,n2,n3) denote 
the fraction of walks which first cross the three barriers after 
(ni, 712,723) steps. Then the excursion set model sets 

pp{p\V) dp = /(M123) dMi23. (22) 

where 

M123 oc (23) 
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denotes the mass associated with the first crossing at 
(721,712,713). We then renormahze the left hand size follow- 
ing the discussion in Section r2.3l In Figure|4]the cyan (solid) 
curves show the predictions of the excursion set approach. 
While they fit better than the perturbation theory's results 
the agreement is not very good at the low density region. 

Although this algorithm allows for different shapes in 
the initial distribution, it only approximates the exact prob- 
lem we wish to solve, because the 7ith step in any given di- 
rection is associated with the same value of Sn, whatever 
the values of n in the other directions. In effect, we are ap- 
proximating 5(711,712,713) with the spherical value S„. for 
each axis i. 

3.2 The ellipsoidal collapse model 

The late time evolution in the ellipsodial collapse model 
is more complicated than in the Zeldovich Approximation 
| |Bond fc Myers 1996[ ). Nevertheless, we can still write the 
nonlinear density as a deterministic (albeit complicated) 
mapping of the three eigenvalues A in the initial field. 

Thus, the non-linear density is given by an expression 
that is analogous to equation (|17p . Namely, 

3 

p=n(i-?o"\ (24) 

1=1 

where in an Einstein-de Sitter universe the evolution of is 
described by 

with 

b. = + and (26) 

j 

This expression for Lext is known as the linear tide approxi- 
mation ( |Bond fc Myers 1996[ ) . The initial condition of equa- 
tion (|25[) is the Zeldovich Approximation. 

The solution of the differential equation above can be 
written as a series: = (^^''^a-' . In what follows we only 
consider the first two terms in this series: 

Cf' = A, (28) 

= ^(/?-2/2) + ^/2 + ^A.(2/i-5A0, (29) 

where Ji = \j and I2 = "l^jy^k ^j'^k (Ohta et al. 2004 
provide a similar expansion). 

We then apply the same argument as in the previous 
section for the perturbation-theory approach to calculate 
the pdf of the evolved density field. Namely, we set 

ppip\V)dp = I d\p{X\l)5n (^P = n(l-'^^')"') ' (30) 

where a^^i = cr^|^' -I- f^^^^ and a includes the shape effect 
(equation [21] with values replacing values). 

Figure [S] shows the correlation between the initial and 
final densities in this model for 8h~^ Mpc cells. Recall that, 
because the evolution depends on all three Ai, whereas the 




5 10 15 20 

P'ni 



Figure 5. Stochasticity in the mapping between linear and non- 
linear overdensity on scales of 8h~^ Mpc. In the upper panel, 
points show the ellipsoidal collapse mapping associated with equa- 
tion l|30| l: solid line shows the spherical collapse mapping of 
equation I I13I I. Both are normalised as described in Section 12.31 
A'^EC = 1.087 for ellipsoidal collapse and A^sc = 1-33 for the 
spherical model. The lower panel shows the difference between 
these two mappings. 



initial density (5l is simply proportional to J]]. A;, we expect 
there to be some stochasticity in p at each (5l. The Fig- 
ure shows that this is not a very large effect. Moreover, the 
solid line shows that the spherical collapse mapping (equa- 
tion [13]) provides a reasonably good description of the mean 
mapping. We exploit this fact in the next section. 

The dot-dashed lines in Figure [4] show the predicted 
pdfs from this second order ellipsoidal collapse model. Note 
that inclusion of these second order terms provides a clear 
improvement over the Zeldovich Approximation prediction; 
this is true over all ln(p) > — 1 in both the 8 and 16h~^ 
Mpc cells. Although the discrepancies at small p are the 
most dramatic (we discuss these in the next section) , we are 
actually more interested in the region around the peak of the 
pdf: the bottom panels show clearly that the Zeldovich Ap- 
proximation curves lie snake around the simulation results, 
whereas this sideways S-shaped residual is largely absent in 
the second order ellipsoidal collapse model. 



4 DISCUSSION 

We discussed two approaches (perturbation theory and ex- 
cursion set) for calculating the dark matter pdf using two 
different models for the dynamics (spherical collapse and the 
Zeldovich Approximation) in each case. Although both ap- 
proaches are deterministic, in the sense that the nonlinear 
evolution is determined by the initial conditions locally, the 
excursion set is slightly less 'local'. We showed that both are 
expected to give similar predictions whenever the dynam- 
ics is approximated by local deterministic mappings (equa- 
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tion[6]), of which the spherical evolution model (equation [13]) 
is a special case (Figure[l]). The agreement is important, be- 
cause the excursion set calculation allows one to connect 
discussions of the dark matter halo distribution with discus- 
sions of the pdf USheth 1998|) . 

Both the perturbation and excursion set approaches, 
when combined with the spherical evolution model, provide 
good descriptions of the nonlinear pdf seen in simulations 
(Figures [2] and |3]). This agreement is slightly unexpected, 
in the sense that the spherical evolution model ignores the 
fact that the nonlinear evolution of a region may be deter- 
mined by quantities other than its initial density. The ellip- 
soidal collapse model is a specific example of this; we studied 
its first (Zeldovich Approximation) and second order expan- 
sions in some detail. In this model, the evolution of a patch 
is determined by its overdensity as well as the surrounding 
shear field (only the overdensity matters for the spherical 
collapse model). As a result, nonlinear evolution changes 
the shape as well as the volume of a patch (in the spherical 
collapse model, only the volume changes). We found that 
including the evolution of the shapes is important when the 
density is high, so this is particularly important when the 
cell size is small. This shape dependence also complicates 
the excursion set calculation (Section 13. ip . 

Despite its increased complexity compared to the spher- 
ical evolution model, we found that the Zeldovich Approx- 
imation resulted in significantly worse agreement with the 
simulations (Figure |4|. Going to second order in the ellip- 
soidal collapse model resulted in more accurate predictions, 
except in the lowest density regions. The relatively large dis- 
crepancies in underdense regions, in both models, may be re- 
lated to the following fact. A Taylor series in 5 will converge 
more rapidly when S > than when S < 0, since the latter 
will be an oscillating series. Thus, it may be that, to ac- 
curately represent the evolution of underdense regions, one 
must go beyond second order in the dynamics. Presumably, 
going to even higher order would further improve agreement 
with the simulations; this is the subject of work in progress. 

Although our perturbation theory calculation accounts 
for the evolution of shapes, our excursion set calcula- 
tion does not. This too is the subject of on-going work. 
The excursion set calculation is particularly interesting 
in view of the fact that the marriage of ellipsoidal col- 
lapse with the excursion set approach appears to provide 
a substantially improved prediction of dark halo abun- 
dances | |Sheth, Mo fc Tormen 200"T] ). Finally, we note that 
our study of the real space pdf can be extended to red- 
shift space; we are in the process of extending previ- 
ous work on this problem (Protogcro s fc Scherrer 1997| 
|Hui, Kofman fc Shandarin 20001 ,Qlita et ITqOO^ '. 

4.1 An extension 

The accuracy of the spherical collapse based predictions, 
and the fact that the associated deterministic mapping from 
initial to final overdensity also provides a good description 
of the mean mapping in the ellipsoidal collapse model (Fig- 
ure[5)l , suggests that one might be able to estimate the shape 
of the initial pdf from a measurement of the evolved one 
as follows. Starting from equation [13] and the normalization 
method in section [231 there is a mapping from the final den- 
sity to its initial value. For each measured M = /oV(1-|-5nl) 



set 



[(l + gNL)/iVsc]~^/'° 

aL(A//Af,c)/<5c 



(31) 



where 5c is the critical linear density associated with col- 
lapse in the spherical model, and A'^sc is the corresponding 
normalization factor (Section 12. 3p . Then make a histogram 
of v, but rather than having each cell count equally, weight 
each by its value of (1 -I- 5NL)/A'ac. The resulting distribu- 
tion of u should provide a good estimate of the shape of the 
initial pdf of (5l/o"l). 

Figure [6] shows the results of this algorithm for a num- 
ber of different scales. Reconstruction of the initial pdf from 
the nonlinear pdf measured on scales larger than 8/i~^Mpc 
works rather well (top panels). The reconstructed distribu- 
tions trace the same Gaussian shape very well; for reference, 
the smooth black curve shows a Gaussian with zero mean 
and unit variance. This mapping works well even when the 
cell size is as small as 4/i~^Mpc (bottom left), but it fails 
at l/i~^Mpc (bottom right). The Figure also shows that it 
is important to account for the factor of 1 -I- 5nl which re- 
lates Eulerian and Lagrangian statistics: weighting each cell 
in the nonlinear distribution equally would have lead one to 
conclude incorrectly that the reconstruction is biased even 
on very large scales. Studying whether or not this is relevant 
for the MAK reconstruction method ( [Mohayaee et al. 2006[ ) 
is the subject of work in progress. 

In effect, the spherical evolution model, when combined 
with the nonlinear pdf, is able to provide a good descrip- 
tion of the initial pdf for scales where the rms ffuctuation is 
smaller than about 2. In principle, this might be used as the 
basis for a test of the Gaussianity of the initial conditions 
because the algorithm makes no explicit assumption about 
the form of the distribution of v. The fact that the non- 
linear pdfs from a range of different scales all map back to 
the same zero mean, unit variance Gaussian curve suggests 
that this method is a good test of Gaussianity. Although 
this is not the first method to reconstruct the shape of the 
initial one-point density field, it is simple and actually fares 
rather well compared to a number of previous methods (e.g. 
Kofman et al. 1994; Narayanan fc Croft 1999). 

This procedure is not as good a test of non-Gaussianity, 
in the sense that it does not yield a simple estimate of the 
shape of the initial distribution if it were non- Gaussian. 
This is because the reconstruction requires an estimate of 
the normalization factor A'^sc, and to calculate it, one must 
first specify the form of the initial distribution. So the al- 
gorithm above is really a self-consistency check: it checks if 
the assumed form of the initial pdf is consistent with that 
reconstructed from the measured nonlinear pdf, under the 
cissumption that spherical evolution is a good model what- 
ever the initial fiuctuation field. Of course, rescaling to u 
involved a rescaling and transformation of the measured Ai, 
but it also involved a{M); this is a consequence of the fact 
that, for Gaussian fiuctuations, information about the vol- 
ume enters only through the scale dependence of the rms 
value. For more general distributions this may not remain 
true, in which case the mapping to the initial variable 5l 
may be more complicated than equation (|3ip . Nevertheless, 
Figure [S] suggests that this may be an interesting avenue 
to explore in future work. We are in the process of extend- 
ing this method to include redshift space effects, as well as 
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Figure 6. Comparison of the reconstructed linear pdf of S/a (histogram) with the expected zero-mean unit-rms Gaussian form (solid 
line). Different panels show results for reconstructions from different smoothing scales. The method works well for large cells (top panels), 
but becomes increasingly inaccurate for smaller cells (bottom panels). When the nonlinear rms fluctuation is larger than 2 (bottom right), 
then the reconstructed pdf is rather different from the Gaussian form, as one might expect from Figure [3] Dotted curves show the bias 
which results from ignoring the extra weighting factor associated with the transformation from Eulerian to Lagrangian statistics. 



galaxy bias. We believe it may be interesting to merge it with 
the method recently proposed by Eisenstein et al. (2007) for 
reconstructing the Baryon Acoustic Oscillation feature in 
the galaxy distribution. 
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where 



In this context, it is interesting that the family of models 
studied by Sheth (1998) included modified Bessel functions 
of the third, rather than second, kind. The solid and dotted 
curves in Figure [1] show equations (|A5|I and (lASP for a few 
values of av- 



APPENDIX A: FULLY ANALYTIC PDF IN 
THE SPHERICAL MODEL 

For a power-law power spectrum (-P(fc) oc fc"), the spheri- 
cal model allows a semi-analytical form of the pdf for some 
special choices of n. To see this clearly, define Sl by 



Sl{p) 



1 



-Sep 



(n + 3)/6 



(1 



-1/Sa 



where av is independent of p. If we set 

Ml' 



av 



6 



then 
P = 

N = 



2<Sc 



2'Kav 



25c 



and the normalised perturbation theory pdf is 



2 , N i^{p/N,5^^ 
P Pip) = „ /„ ^ exp 



6 



where 

K.{p,5c)=p' " + p ' " 

Setting 5c — 5/3 means n — —6/5 and so 
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(Al) 
(A2) 

(A3) 
(A4) 

(A5) 
(A6) 
(A7) 



were K„{x) is a modified Bessel function of the second kind. 
The random walk prediction for this model is 



P P{p) 



(P/N) 



3/10 
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